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In this paper intermodal modulational instability (IM-MI) is analyzed in a multimode fiber where 
several spatial and polarization modes propagate. The coupled nonlinear Schrodinger equations 
describing the modal evolution in the fiber are linearized and reduced to an eigenvalue problem. As a 
result, the amplification of each mode can be described by means of the eigenvalues and eigenvectors 
of a matrix that stores the information about the dispersion properties of the modes and the modal 
power distribution of the pump. Some useful analytical formulas are also provided that estimate 
the modal amplification as function of the system parameters. Finally, the impact of third-order 
dispersion and of absorbtion losses is evaluated, which reveals some surprising phenomena into 
the IM-MI dynamics. These outcomes generalize previous studies on bimodal-MI, related to the 
interaction between 2 spatial or polarization modes, to the most general case of IV > 2 interacting 
modes. Moreover, they pave the way towards the realization of wideband multimode parametric 
amplifiers. 


I. INTRODUCTION 


The past 40 years have witnessed a huge increase of the 
transmission capacity in single-mode fibers. Every avail¬ 
able degree of freedom has been explored, such as mul¬ 
tiplexing in time, in wavelength, in polarization and in 
phase. On the other hand, nowadays single-mode fibers 
are gradually reaching their capacity limit of about 100 
Tbit/s which is dictated by the Shannon theorem |2j. 
In this scenario, multicore and multimode fibers are 
rapidly emerging as the ideal solution in order to fulfil 
the growing demand of capacity. In these fibers nonlinear 
pulse manipulation can be efficiently exploited for various 
applications [3]. Furthermore, they allow exploring a fur¬ 
ther degree of freedom, namely the spatial, and therefore 
the implementation of space division multiplexing (SDM) 
schemes [l] . In SDM each fiber core and/or mode repre¬ 
sents an independent information channel parallel to the 
others, which permits in principle to largely overcome the 
capacity limit of single-mode fibers. 

A key-component for an efficient operation of SDM 
schemes is represented by multimode optical amplifiers 
[5j, which have been the subject of intense research over 
the last years. Among the most important we find multi- 
mode fiber Raman amplifiers ( RA) [Tj and multimode 
erbium doped fiber amplifiers (EDFA) [8/ 0]. 

An alternative fiber-based approach is represented by 
parametric amplification m, which is based on a mod¬ 
ulational instability (MI) process. 

Compared to RAs and EDFAs, parametric amplification 
guarantees more flexibility since it makes possible to tune 
the position of the amplified sidebands at arbitrary wave¬ 
length by means of the input pump and of the fiber pa¬ 
rameters. Moreover, Kerr-nonlinearity is practically in¬ 
stantaneous, enabling ultrafast signal processing applica¬ 
tions. 

Intermodal-MI (IM-MI) in multimode fibers was demon¬ 
strated decades ago m and more recently explored as 
a possible means for supercontinuum generation m as 


well as for the development of efficient fiber parametric 
amplifiers [HE] and low-cost broadband fiber sources 

Ca¬ 
in all previous works, IM-MI has been analyzed by con¬ 
sidering separately the interaction between couples of 
modes and studying the bimodal-MI process related to 
the phase-matching of each couple. With this approach, 
the complex IM-MI dynamics is decomposed in a set of 
bimodal-MI processes, which greatly simplifies its analy¬ 
sis. On the other hand, the presence and the influence of 
the non-phase-matched modes is usually neglected. 

The main purpose of this paper is to go beyond this sim¬ 
plistic view and to develop a model taking into account 
for the interaction between all the fiber modes, which is 
the solely way to correctly describe the IM-MI. Further¬ 
more we provide, for the first time to the best of our 
knowledge, some useful analytical estimates for the IM- 
MI gain and we evaluate the impact of higher-order dis¬ 
persion terms and losses, which reveals some surprising 
phenomena into the IM-MI dynamics. The final result is 
a general description of the IM-MI that paves the way to 
the realization of wideband multimode parametric am¬ 
plifiers. 

The paper is organized as follows. 

In Section [TT] we linearize and reduce to an eigenvalue 
problem the coupled nonlinear Schrodinger equations 
(CNLSE) describing the propagation in a multimode 
fiber. We focus our attention on single-core isotropic, 
highly-birefringent and telecommunication fibers. We fi¬ 
nally obtain a matrix M storing all the information about 
the dispersion properties of the modes and the modal 
power distribution of the pump. This outcome general¬ 
izes previous results on bimodal-MI [TB] , generated by 
the interaction between 2 spatial or polarization modes, 
to the most general case of N > 2 interacting modes. 

In Section |III| we find some useful analytical formulas 
that quantify the growth of each mode as function of the 
eigenvalues of M with negative imaginary part, namely 
the gains, and their corresponding eigenvectors. Inter- 
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esting enough, several gains could coexist for the same 
pump-sideband detuning, giving rise to a competition be¬ 
tween two or more amplification processes. We thus in¬ 
troduce the concept of dominant gain, which is the truly 
leader of the MI growth, and we show that different gains 
could be dominant at different fiber positions and for dif¬ 
ferent modes. 

In Section |IV| we provide a deeper physical insight into 
the IM-MI and we derive some analytical estimates for 
the IM-MI gain as function of the pump-sideband detun¬ 
ing and the system parameters. We put in evidence that 
IM-MI is truly associated to a phase-matched four-wave 
mixing process, which is the usual point of view of past 
works, but we also clearly show that all the fiber modes 
undergo amplification, and not only the phase-matched 
ones. The amplification of each mode can be controlled 
by means of the system parameters and over a large band¬ 
width that increases with the number of the modes. 

In Sections [Vj [VI] and US we analyze the impact of 
third-order dispersion as well as of the pump power dis¬ 
tribution and of the propagation losses on the IM-MI 
dynamics. We find that absorption losses cause a shift 
of the gain curves towards low frequencies and that, dif¬ 
ferently from single-mode fibers, third-order dispersion 
plays a substantial role and may lead to the generation 
of secondary MI bands in addition to conventional ones. 
In Section |VIII| we discuss the limits of validity of our 
model and finally, in Section |IX[ we resume the main 
outcomes and we give the concluding remarks. 


II. GENERAL THEORY 

Let us consider a beam, centered at the carrier angular 
frequency ui p , which is injected in an optical fiber and is 
coupled to the N propagating modes of the fiber. In the 
weakly guiding approximation, which is well accurate in 
typical silica fibers with low core-cladding refractive in¬ 
dex difference, the modal transverse profiles M n (x,y) ( 
1 < n < N) are linearly polarized and almost indepen¬ 
dent of their polarization. The a—polarized component 
(a = {x, y}) of the total electric field E propagating in 
the fiber along the z-axis can be written as follows: 


N 

E a = ^2 Pna{z, t)exp(i/3 p , na z + iui p t)M n (x , y) + c.c( 1) 

n—1 

where P na is the envelope of the a-polarized n-mode, 
which will be called na —mode throughout the paper, and 
Pp,na = Pna(u p ) is the corresponding propagation con¬ 
stant at the frequency lo p . 

Starting from the Maxwell equations, with a polarization 
that takes into account the nonlinear cubic response of 
silica, the following set of coupled equations can be de¬ 
rived which describes the spatio-temporal dynamics of 
the a:-polarized modal envelopes mi: 


dP nx _ 

dz 

— l &Pnx . $2nx & Pnx fisnx & Pnx p 

~ Vnx ~dT ~ + ~6 di~ a ™ 

i ^ ^ C-klmn (PkxPlxPmx^^-,klmn 4 “ ^P^ x Pi x Pmx^2,klmn) 4 “ 

klm 

i ^ ^ Cklmn i^PkyPlyPmx^'^iklmn 4 “ ^ PkyPlyPmx^4,klmn) 

klm 


( 2 ) 


An equation similar to Eq. § is valid for the y- 
polarized modal envelopes P ny by exchanging the labels 
a ' y. 

Here * denotes the complex conjugate and e r ^imn = 
exp(iA/3 r ,kimnZ ), r = {1,2,3, 4}, being: A^i Mrnn = 
fip,kx 4“ Pp,lx -P P ,mx -P P ,nx 5 A 02 ,klmn — ftp ,kx 4“ ftp,lx 4“ 
Pp ,mx - P P ,nxi A/? 3 ,klmn — Pp,ky T Pp,ly Pp 7 mx Pp^nxi 
^P^^klmn — Pp,ky T Pp,ly T Pp,mx Pp,nx' 

Coefficients v na = du/dp na \ Up , /3 2na = d 2 /3 na /du 2 |^ p , 
p3na = d 3 P na /du 3 \ bJp and a na indicate respectively the 
group velocity, the group velocity dispersion (GVD), the 
third-order dispersion (TOD) and the propagation losses 
related to the na-mode at the frequency lo p . 

The nonlinear coupling coefficients read as Ckimn = 
( n 2 uj p /c ) J f MkMiM m M n dxdy , where n 2 is the nonlin¬ 
ear index of the fiber, c is the speed of light in vacuum 
and all the modal profiles are normalized so that the area 
./' I M 2 dxdy = 1. 

Equation (|2j) is valid in a fiber where we can neglect 
the linear coupling among the different modes, which is 
induced by random perturbations in the fiber structure 
such as manufacturing imperfections, environmental vari¬ 
ations or local mechanism stress PT5j ■ 

For example, in a highly birefringent (HiBi) fiber the 
intrinsic large birefringence makes the propagation con¬ 
stants of all spatial and polarization modes to be notice¬ 
ably different so as to nullify in practice the linear modal 
coupling. 

On the contrary, in quasi-perfectly circular core fibers 
we find group of modes that are quasi-degenerate and 
may thus experience a strong linear coupling. Neverthe¬ 
less Eq.Q is still valid provided that the length L of 
the fiber is small if compared to the characteristic lin¬ 
ear coupling length of these modes. Actually, under this 
condition the fiber may be considered isotropic so that 
quasi-degenerate modes do not exchange energy and can 
thus be represented by an unique mode that is a proper 
linear combination of them. 

We introduce here the nonlinear length L^l = (yT) -1 , 
where 7 = max{ckimn} is the largest among the nonlin¬ 
ear coupling coefficients and T = I Pnx | 2 + I-Pny 1 2 is 
the total input power. This parameter provides an idea 
of the length scale of the nonlinear interactions and per¬ 
mits to select in Eq. § the only relevant nonlinear terms, 
that are those for which the condition \&.p r ,kim n LNL\ ~ 0 
applies. Otherwise, if \A/3 rt kimnLNL\ » 0, then the 
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corresponding nonlinear term is rapidly oscillating and 
averages out to 0, so that it can be neglected. A proper 
discussion about the validity of this rule of thumb is done 
in Sec. IVIIII 

In practice, due to the large phase-mismatch between the 
different spatial modes, in a single-core HiBi or isotropic 
fiber the condition \Af3\ tk imnpNL\ ~ 0 is achieved only 
if {k = mj = n} or {k = n,l = to}. Similarly 
{k = to, l = n} or {k = l, m = n} are needed in or¬ 
der to have |A/? 2 ,kimnpNLl ~ 0. In a HiBi fiber, where a 
large phase-mismatch is present even between two polar¬ 
ization modes, the condition \Ap^^ k i mn LxL\ ~ 0 cannot 
typically be attained, while \Ap^ k i m nLNL\ ~ 0 requires 
{l = fc, to = n}. 

Taking into account only the relevant nonlinear terms, 
Eq .([2]) can be rewritten as follows in the case of an HiBi 
fiber: 

dPnx _ -!dP nx .f3 2 nxd 2 P nx fcnx d 3 P nx 

dz ~ Vnx dt 1 2 dt z + 6 dt 3 

^■nx-^nx + 'ibsC'nnl Pnx | Pnx “1“ ^ ^ ^ ^ | | C kn | Pkx | Pnx 

k^n 

+ ibxC n n\Pny\ Pnx *E ^A-P'kn\Pky\ Pnx (3) 

k^n 

where C kn = 3 c kknn ; = 1 is the coefficient related 
to the self-phase modulation (spm); bx = 2/3 is the coef¬ 
ficient related to the intramodal cross-phase modulation 
(xpm); 6 || = 2 and b± = 2/3 are the coefficients related 
to the intermodal xpm involving modes with parallel or 
orthogonal polarization, respectively. 

Equation Q holds also for an isotropic fiber in the 
scalar case, that is when the input field is linearly po¬ 
larized along one fixed direction, let us say the a:-axis 
(i.e. we set bx = b± = 0). If the input field at 
the isotropic fiber is right or left circularly polarized, it 
proves convenient to rewrite Eq. © in terms of the cir¬ 
cularly polarized components Pn+ = (Pnx + iPny)/V 2 
and P ra _ = ( P nx — iP ny )/V 2 , which brings to the same 
equality of Eq.Q after substitution of label x (y) with 
+ (—) and setting bs = 2/3, 6 || = 4/3, bx = b± = 0. 
Differently from HiBi or short isotropic fibers, in telecom¬ 
munication fibers rapid random variations of the fiber 
structure must be taken into account, as they lead to a 
non-negligible linear modal coupling. At this purpose, a 
generalized multimode Manakov model has been derived 
by Mumtaz et al. in m- It turns out that in the CW 
limit Eq. Q is still valid in a telecommunication fiber by 
setting bs = bx = 8/9, b\\ = b± = 4/3 and replacing 
PnxjPny by the envelopes P nx ,P n y that are obtained by 
means of a proper unitary transformation. 

In the following we will indicate the field components as 
Pnx ( P n y ), with the assumption that they should read as 
P n+ ( P n _) in the case of isotropic fibers where right/left 
circular polarization are involded, and as P nx ( P ny ) in 
the case of telecommunication fibers. 

As usual, in order to examine modulational instability 


processes, we introduce small amplitude perturbations. 
We thus decompose the modal envelope P na (a = {x, y}) 
in the sum of a pump p and two perturbations s and i 
that indicate a signal Stokes and an idler anti-Stokes side¬ 
band symmetrically detuned with respect to the pump 
and centered at the frequencies uj s and w*, respectively: 

Pna = Pna + S na eXp(iAP$>*) z + lflt) + 

i na exp(iA/3f 3> f' ) z - iCtt) (4) 

where H = ui s — u> p = co p — uii is the pump-sidebands 
detuning. The wavevector mismatches read as A/%^ = 

fis,na [3p,na £md A {3na^ = f^i,na Pp,na> where fis,na = 

Pnai^s) and pi^na = /3 na {ui) indicate the propagation 
constant of the na -mode at the frequencies lo s and u>i, 
respectively. 

The decomposition Eq.Q is inserted in Eq.Q in order 
to find the set of differential equations that rules the evo¬ 
lution of pump and sidebands. 

The equation for p nx is given by Eq. © after substitu¬ 
tion of P with p. We initially neglect propagation losses, 
namely we set a nx = a ny = 0. The system dynamics 
in the presence of losses will be discussed later (see Ap¬ 
pendix 2). Under this condition and in the CW-limit one 
obtains the following analytical solution: 

Pnx(z) = \pnx\exp(iip nx (0) + i(j> nx z) 

4*nx ^sP'nn | P nx | H - bxCnn\Pny\ “I” 

+ E kiicwbfczh + ^2 bj_C k n\Pky\ 2 (5) 

k^n k^n 

where is the phase of the input p nx { 0). A simi¬ 

lar analytical solution holds true for p ny after exchanging 
x -H- y. 

The equation for the modal envelope s nx , obtained by 
linearization of Eq.Q and neglecting losses, reads as: 

——— — ii^(f)nx T bsCnn \Pnx | )Sna;~f“ 

OZ 

ibxC nn p* n yP nx s ny exp{iA^ s) z - iA/3^ s) z)+ 
ib s C nn pl x i* nx exp(-iA/3^ s) z - iA^ l) z)+ 
ib x C nn p nx p ny il y exp(-iA^ s) z - iAf3 { n P3) z)+ 
i ^2 b\\C kn pl x p nx s kx exp{iAl3 { ^ s) z - iA/3^ s) z)+ 

k^n 

i ^2 b±C k nPkyPn X Sk y exp(iA/3^y S) z - iA/3^ s) z)+ 

k^n 

i ^2 b \\ c knPkxPnxi* kx exp{-iA^' l) z - iAf3^ s) z)+ 
k^n 

i E b±C knPkyPnxil y exp{-iA^ P 2 )z ~ z) (6) 

k^n 

Similar equations apply to s ny after exchanging x -H- y , 
to i nx after exchanging i o s and to i ny after exchanging 
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x -o- y, i -o- 5 . 

In order to get rid of the oscillating terms, it proves useful 
to employ the following change of variables: 


Sna(z) = S na (z)exp(ilp na ( 0 ) + i4>naZ - lA/3^z) 
ina(z) = i na (z)exp(ilp na ( 0 ) + i(f>naZ ~ iA/3^ z) 

(7) 

After insertion of Eq. 0 in Eq. 0 and rewriting the 
fields Pnx,Pny according to the solution given in Eq.([5]), 
we finally obtain the following eigenvalue problem: 


cbv = iMv 


v = [s x s y i x * i y *] T 


-^Isx.sx 

M SX:Sy 

M sxix 

M SXjiy ' 

-M-sy,sx 

Mgy.sy 

M sy . ix 

M S y,iy 

-M ix?sx 

—M ixsy 

-M ix! ; x 

-Mix.iy 

-M iy , sx 

-Miy.sy 

-M iy>ix 

-Miy.iyj 


( 8 ) 


where s a = [s ia s 2a •■•SiVaF and i a = [i la i 2a —i N a} T 
(a = {x,y}) are N x 1 vectors, and M is a 4iV x 41V 
matrix composed by 16 N x N matrix blocks. 

Each matrix block accounts for the interaction between 
different subsets of sideband modes. For example, the 
block M sx sx is related to the mutual interaction between 
the x —polarized signal modes, whereas M sx sy is related 
to the interaction between the x— and y —polarized signal 
modes. The elements of M sx sx , M SX!Sy and M SX! i x are: 


M 

M 

M 

M 


SX,SX 

SX,SX 

sx,sy 

sx,sy 


M s 

M s 


[n,n] = A/3^ s) + b s C nn \p nx \ 2 
[n, m] = b\\ C mn | p mx 11 Pnx I ( n^m ) 

[n, 7T-] — bxCnn \Pny \ \Pnx \ 

[n,m\ = bj_Cmn\p m y\\p n x\ (n ^ m) 
[ 77 ., 7z] = bgC nn \Pnx \ 

[n, m] = 6|| C mn I Pmx 11 Pnx I (n^m) 


(9) 


Furthermore Al sx i y — fM sy i x — JVL , x .i y — Nl sxsy , 
the matrix M syiSy reads as M sx ,s X after replacing x with 
y, M ix , ix reads as Mg XjSX after replacing s with i: Mi y i y 
reads as M syj s y after replacing s with i: M sy i y reads 
as Ms X .i x after replacing x with y ; M sy:SX = M SXi g y T ; 
M ix , sx = M sx , ix ; M ix , sy = M syiix T ; M sy , sx = M SX!Sy T ; 
M iy , sx = M SXj i y T ; M iy , sy = M sy , iy ; M iy , ix = M ix , iy T . 
The solution of Eq.Q reads as: 


4 N 

= E c k w k [j]exp(i\ k z) (10) 

fe =1 

where X k and w k (1 < k < AN) are respectively 
the eigenvalues and the eigenvectors of M, whereas v\j\ 


(■ w k [j ]) indicates the j —element of v (wj,). The coeffi¬ 
cient Cfc is given by the projection of the input v(z = 0), 
which is fixed by the input modal envelopes of the side¬ 
bands, over the set of eigenvectors. We highlight that the 
structure of the matrix M guarantees that its eigenvec¬ 
tors satisfy the orthogonality relation w*,i • Dwt 2 = 0 
(ki ^ k 2 ), where D = diag([ones 2 N, —ones 2 N ]) is a di¬ 
agonal matrix, ones 2 N is a 1 x 2N vector of ones and 
• indicates the scalar product. By exploiting this or¬ 
thogonality relation we can easily compute c k as follows: 
Cfe = (v(2 = 0) • Dw k )/(w k • Dw k ). 

An eigenvalue X k with negative imaginary part leads to 
an exponential amplification of v[j] which is character¬ 
ized by a gain g k = — Im(X k ). We may be tempted to ap¬ 
proximate Eq.( 10) by taking into account only the largest 
gain g max and its corresponding eigenvector w max , as 
they ultimately determine the asymptotic evolution of 
v\j], namely \v\j](z)\ « \c max \\w m a X \j]\exp{g m axz). 
Nevertheless, in doing so, we may not properly de¬ 
scribe the evolution of v[j] when two or more eigenval¬ 
ues with negative imaginary part are present. Indeed 
two gains gi and g 2 could exist such that g 2 > g\ but 
\ciWi[j)\exp(giz) > \c 2 w 2 [j]\exp(g 2 z) at a certain posi¬ 
tion z. In this case we will define g\ as the dominant 
gain for v[j] at the position 5, as it plays the role of true 
leading term in the MI growth of v\j\. 

From what discussed above, a good approximation of 
Hj]| at the fiber exit, let us say in z=L, could be ob¬ 
tained by computing the corresponding dominant gain, 
that is \v\j]\(L) « max k {\c k w k \j]\exp(g k L)}. The signal 
and idler modal amplitudes at the fiber exit can be finally 
estimated by noting that |s na; | = |t>[n]|; |s ray | = | 7 i[n+A^] |; 

| inx | = | v[n + 27V]|; \i ny \ = \v[n + 37V] | (1 < n < N). 

We underline that the coupling coefficients C kn of 
the matrix M are generally fl-dependent, and so are 
the wavevector mismatches A/Fa’ and A/Fa’F Both 
the mismatches could be approximated by expanding 
/3na(<x s ) and f3 n a(wi) in a Taylor series centered around 
Pna(iXp)- Truncating the expansion at the third-order we 
finally get: 


A/3 nci S ' > — ^/Vna + Xl 2 p 2na /2 + f2 3 p3na,/6 

A = -n/v na + fi 2 /W2 - n 3 p 3na /6 ( 11 ) 

Therefore for a given detuning and knowing the 
input modal envelopes of both pump and sidebands 
we may compute the corresponding eigenvalues and 
eigenvectors of M and then estimate the output modal 
envelopes s na (L) and i na {L) as formerly discussed. 

From what stated it clearly appears that the MI growth 
of the sidebands depends on the eigenvectors and 
eigenvalues of M, therefore it can be controlled by 
means of the system parameters, such as the modal 
power distribution of the pump, the nonlinear coupling 
coefficients and the modal dispersion characteristics 
(that is group velocity, GVD and TOD of each mode). 
This outcome opens the way towards parametric amplifi¬ 
cation in multimode fibers and represents the core of the 
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IM-MI analysis developed in this paper. It generalizes 
previous results on vectorial MI in single-mode fibers, 
for which the matrix M is 4 x 4, to the most general case 
of the IM-MI in a fiber where N spatial and polarization 
modes propagate. Furthermore, as previously observed, 
the theory developed in this section applies to all the 
main types of commonly used optical fibers. 


III. AMPLIFICATION OF AN INPUT ASE 
NOISE 

In order to provide numerical evidence confirming the 
theoretical results previously argued, we simulate the 
propagation in a multimode fiber by numerically solv¬ 
ing Eq. © with the split-step Fourier method. Finally, 
we compare the output modal amplitudes obtained from 
numerical simulations with those inferred from the ana¬ 
lytical outcomes discussed in Section [n] 

With the aim of focusing on the most striking features of 
the IM-MI, we represent a simple but effective example: 
the propagation in an isotropic fiber of 4 spatial modes 
that are linearly polarized along the same direction. 

The simulated fiber is circular in cross-section and 
step-index, which allows easily solving the well-known 
modal characteristic equation in order to compute the 
modal transverse profiles and the main modal parame¬ 
ters (Fig. [l]and Table 1-2). 

The CW input field is centered at the carrier wave¬ 
length A = 1550 nm and is linearly polarized along the 
a;—direction. The diameter of the core is 2 R = 24/xm, 
the refractive indexes of cladding and core are n c i a d = 1.5 
and n core = 1.5035, respectively. The corresponding V- 
number of the fiber V = 2nRX~ 1 (nl ore — n ): ia[i ) 1//2 is 
V ~ 5, so that 4 non-degenerate spatial modes propa¬ 
gates : LP 0 i , LP 0 2 , LPu and LP 2 \. In our case these 
modes are I—polarized and will be indicated respectively 
as lx—mode, 2x—mode, 3x—mode and 4x—mode. 

We don’t take into account for the 2-fold degeneracy of 
the LPu mode: indeed we assume that only one among 
the two degenerate modes is excited at the input fiber 
and that it does not couple to the second because the 
fiber is ideally isotropic. A similar consideration applies 
to the LP 2 i mode. 

We set the total input power equal to 4000IU and 
uniformly distributed over the 4 pump modes, there¬ 
fore \p nx \ 2 = 1000W (n={l,2,3,4|). Such an high- 
power injection is indeed feasible in multimode 
fibers when employing nanosecond pump sources 
at low rates, so that the pump can be practically 
considered CW [I9l '20]. We also add a weak back¬ 
ground white noise, which could be the amplified spon¬ 
taneous emission (ASE) noise in a realistic experiment. 
The total input field for the nx— mode is thus written as 
P n x(0,t) = Pnx (0) +r n x (0, t), being p nx ( 0 ) = ( 1000 W) 1/2 
and r nx (0,t) a white noise generated by adding at each 
frequency component a random variable with indepen- 
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FIG. 1. Modal transverse profile of the 4 modes supported by 
the simulated fiber. The radius of the core is R = 12p.m. 


TABLE I. Propagation constant /3, group velocity mismatch 
( GVM ), group velocity dispersion j3 2 and 3rd order disper¬ 
sion @3 of the 4 supported modes at the pump wavelength 
1550nm.These parameters have been adjusted in order to 
limit the IM-MI bandwidth to nearly 30 THz, so to avoid 
higher-order dispersion terms besides the TOD. Note that in 
a real fiber they can be adjusted over a wide range of values 
by means of an appropriate fiber design. 



LP 01 

LP 02 

LPu 

LP 21 

P[pm *] 

6.0995 

6.0836 

6.0891 

6.0848 

GVM[psm~ L ] 

0 

10.8 

7.1 

13.6 


21.7 

-147.7 

36.3 

-3.5 


89.5 

-7361.1 

-169.9 

-2128.7 


dent Gaussian-distributed real and imaginary parts. 

We finally solve Eq. © by split-step Fourier method us¬ 
ing P nx (0,t) as input field for the nx—mode. 

Each pair of noise frequency samples centered at +0 and 
—f 1 play the role of small input signal and idler perturba¬ 
tions that are amplified by the IM-MI process. The power 
spectrum R nx (L, 0) of r nx (L,t) reveals the amount of 
power carried by the nx—mode at the detuning O and 
at the fiber exit z = L. As explained in Appendix 1, we 
make use o f the averaged power spectrum R nx avg (L,Q) 
(Fig. Al-1), which expresses the truly noise power level 
around SI. 

The function R n x,avg(L : fl)/R nXtavg (0,fl ) defines the ra¬ 
tio between the output and the input power coupled 
to the nx-mode, and the logarithmic square-rooted 
power-ratio normalized to distance, here indicated with 
A nx (L, f2), represents the amplitude amplification factor 
for the nx—mode at the fiber output: 


Ana:(A, If) — ^Log 


R n 


;(A, fl ) 1/2 ' 


R 


•nx,avg 


( 0 , O) 1 / 2 


( 12 ) 


We will compute the function A nx obtained from nu- 
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TABLE II. Nonlinear coupling coefficients C k n of the sim¬ 
ulated fiber. The coefficients are normalized with respect to 
Cn = 10 W~ 1 km~ 1 , so that the element in row k and column 
n represents Chn/Cu- 



n = 1 

n = 2 

n = 3 

n = 4 

k = 1 

1.00 

0.73 

0.66 

0.45 

k = 2 

0.73 

0.96 

0.37 

0.33 

« 

II 

0.66 

0.37 

1.04 

0.61 

k = 4 

0.45 

0.33 

0.61 

0.92 


nonlinear time. According to this normalization, in the 
example under discussion £ = 1 corresponds to z = 
0.1m and v = 1 corresponds to a detuning (27r) _1 f2 = 
30.IT Hz, whereas a normalized gain g = 1 corresponds 
to a real gain of 10 m -1 . 

In Fig. [2]the gain coefficients are displayed as function of 

v. 

We notice the existence of several gain curves, which is 
due to the rich set of modal interactions. This point will 
be treated in detail in the next section. 


merical solution of Eq.([3]) and we will compare it with 
the following analytical estimate A nx ^ est gathered from 
the results exposed in the previous section (see Appendix 
1): 

-Anx,est (Z/, — 9dom,nx H - L Log{\Wdorn,nx\p\\) 

Wdom.nx = Wdom.nx/(^dom.nx • DWd om ,nx) 

9dom,nx = 9dom,nx (L,n) w do 

m,nx — Wdom.nx (L,n) 
(13) 

where gdom,nx is the dominant gain for the nx—mode 
at the detuning Q and at the position z = L, whereas 
Wdom.nx is its corresponding eigenvector, and they sat¬ 
isfy the following condition: 


\wdom,nx[n]\exp(gdom t nxL) = max k {| Wk [n]\exp{g k L )} 

(14) 


Note that the third line of Eq. (13) has been added 
to emphasize the dependence of gdom,nx and Wdom.nx 
on both L and f2 and that Eq. (13) is valid whether 
A nx ,est{L, A2) > 0 , otherwise we set it equal to 0 . 

In order to find gdom nx and Wdom.nx we build the 16 x 16 
matrix M of Eq.([ 8 | by setting \p nx \ = (10001F) 1 / 2 , 
\Vny\ = 0 (n={l,2,3,4}) and using the modal parameters 
and coupling coefficients displayed in Table 1-2. At each 
detuning f 2 we then compute the corresponding gain co¬ 
efficients g k = —3(A*,) and eigenvectors Wk (1 < k < 16) 
and we look for the couple (g k , Wk) that fulfill Eq. (14). 
We point out that the structure of the matrix M guar¬ 
antees that 3 fc(fl) = g k (—Q). Furthermore the idler com¬ 
ponent w k [n + 21V] is directly tied to the signal compo¬ 
nent w k [n] by the relation w k [n + 21V](12) = w k [n](—fl). 
For this reason we compute only w k [n] as function of 
12 , assuming that positive frequencies refer to the sig¬ 
nal nx—mode, whereas negative frequencies to the idler 
nx—mode. Note also that an equation similar to Eq. (13) 


would apply for the amplification factor A nytest (L, 12), 
here not considered as the modes are x-polarized, after 
replacing of x with y and of w k [n] with w k [n + IV]. 

For convenience in the following we make use of the nor¬ 
malized dimensionless distance £ = z/Lnl ,i and de¬ 
tuning v = (27t)^ 1 12Ttvl,i, where is the non¬ 

linear length related to the lx—mode and Tnli = 
(l/? 2 , lx \Lnl,i/2) 1 ^ 2 is the corresponding characteristic 



FIG. 2. normalized IM-MI gain VS normalized frequency de¬ 
tuning v. Gain values are normalized with respect to the 
nonlinear length. Different gain curves are found which are 
plotted with different colors and labeled with Ba, Bb, Bj. 
Panel (a): global view. Panel (b): zoom in the band from 
v = 0.85 to v = 1.1.Panel (c): zoom in the band from v = 0 
to v = 0.65. Note that the IM-MI gain is symmetric with 
respect to v = 0, therefore the same gain curves are found at 
positive and negative detunings. 


Here we highlight that some of the gain curves overlap, 
that is to say, for some value of v two gains coexist 
that are dominant at different fiber positions. Without 
any loss of generality, we assume that the two gains 
come from the first two eigenvalues of M: therefore 
we indicate with g\ = —S(Ai) and 32 = ^S(A 2 ) the 
concurrent gains and with w 1 and w 2 the corresponding 
eigenvectors . 

According to Eq. (14) and assuming g 2 > gi, if 
|wi[n]| > |tZ> 2 [xz.] | then the gain gi is dominant, for the 
nx—mode, whenever \wi[n]\exp{giz) > \'W 2 [ r >]\^xp(g 2 z), 
that is when z < (h — h)/(g 2 — gi), being 
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l 1 = Log{\w\[n]\) and I 2 = Log(\ui 2 [n]\). Other¬ 
wise, if |u)i[n]| < |"iZ >2 ["Ti] |, then g 2 is dominant at any 
fiber position for the nx— mode. 

As example, let us consider the gain curves Bp and Bq, 
which overlap in the whole band of Bq (see Fig. [ 2 ]). 
Figure [3}i shows a zoom of Bp and Bg around the 
detuning v = —0.43. We observe that for v = —0.43 
the peak g\ = 0.71 of Bq is lower than the gain value 
g 2 = 0.90 of Bp ; nevertheless Log(\u>i{2]\) ss —0.35 is 
larger than Log(\w2[2]\) ~ —3.35 (Fig. [3 o), which makes 
gi to be dominant for the 2x— mode until £ = 15.8. More 
generally we can conclude that for any detuning close to 
v = —0.43 the gain curve Bq, and not Bp, is the most 
appropriate for describing the early-stage growth of the 
2x— mode. 

This is confirmed by the numerical results depicted 
in Fig. |4j where the amplitude amplification fac¬ 
tor A 2x (L,v) obtained by numerical solution of 
Eq. © is displayed when the fiber length L = 5. 
In Fig. ©’ we observe that, around v = —0.43, 

A 2 x(L,v) « Bg{v) — 0.07, that is to say A 2x is a copy 
of Bq lowered of —0.07 = L~ 1 Log{w i[2]), as predicted 
by Eq.(fl3|). 


On the contrary, when L > 15.8 the gain g 2 is dominant 
for the 2a;—mode at v = —0.43: the gain curve Bp 
is therefore the most appropriate for describing the 
second-stage growth (i.e. for £ > 15.8) of the 2a;—mode. 
This is confirmed by Fig. [5] where the amplification 
factor A 2 X {L,u) is shown when L = 16; in this case, 
around v = —0.43, A 2x {L,v) ss Bp(y) — 0.21, with 
—0.21 = L~ 1 Log(w2[2]) as predicted by Eq.(13). 
Differently from the case of the 2a;—mode, the gain 
g 2 is always dominant for the 4a;—mode in the band 
of B g centered at v = —0.43, because in that band 
aJ) 1 [4] < fib [4] (see Fig. §). The numerical solution of 
Eq.® confirms indeed that the gain curve Bp is the 
most suitable at describing the IM-MI growth of the 
4a;—mode even at early-stage ( see Fig. [6]). 


The results displayed in Figs. [4]- [6] show a good agree¬ 
ment between the numerical simulations and the analyti¬ 


cal estimation Eq. (13), which represents therefore a sim¬ 


ple and powerful tool for describing the IM-MI growth of 
the modes. 

In Fig. [4] we see that the 2x—mode is amplified in the 
bands of all the gain curves Ba...j , which confirms what 
previously observed, namely that each mode of the fiber 
is amplified in the whole IM-MI band and its amplifica¬ 
tion depends on the eigenvalues and eigenvectors of M, 
which in turns depend on the system parameters. 

We also point out that the term L -1 Log(\uidom,nx[n]\) 
in Eq. (13) explains the strong asymmetries in the am¬ 


plification spectrum of the fiber modes which can be 
observed comparing panel (b) with (c) or (a) with (d) 
in Fig. [4] In fact the dominant gain at the signal and 
idler detuning is the same, i.e. gdom,nx(L,v) = 

gdom,nx{L, —u), but the corresponding eigenvector com¬ 
ponents are in general different, i.e. Wdom,n X [n]{L,v) ^ 


Gain 



Log( |w 1 2 [2]|) Log( |w 1 2 [4]|) 


0(b) 


Log( |w,[2]|) 


-1 

-2 

-3 

-4 


Log( |w 2 [2]| ) 


-0.43 

v 


0(C) 


Log( |w [4] |) 


c5g( |w [4]| 


—0.43 

V 


FIG. 3. Panel(a): IM-MI gain around v = —0.43. We rec¬ 
ognize the gain curves Bg (red) and Bp (magenta). The 
black dots identifies gi = Bg(u = —0.43) = 0.71 and gi = 
Bp(i/ = —0.43) = 0.9. Panel(b): eigenvector components 
related to the 2x-mode: Log(|iZ>i[2]|) (red) and Log{\u> 2 [2]|) 
(magenta).The black dots identifies Log(\w\[2]\) = —0.35 and 
Log( |iZ >2 [2] |) = —3.35 computed at the detuning u = —0.43. 
Panel(c): eigenvector components related to the 4x-mode: 
Z/Ogr(|uii[4]|) (red) and Z/Ogr(|tI> 2 [4]|) (magenta). The black 
dots identifies Lo(;(|«)i[4]|) = —3.22 and Log(|W 2 [4]|) = —0.35 
computed at the detuning v = —0.43. 


vJdom,nx[n](L,-v ), which makes A nx {L , v) ^ A nx (L, -v) 
according to Eq.(13). 

To resume: the key concepts that should be retained from 
this section are essentially two. First, the modal amplifi¬ 
cation may exhibit a complex dynamics, which is due to 
the presence of different gains that are dominant at dif¬ 
ferent fiber positions. Second, the amplification of each 
mode, well described by Eq. (13), can be controlled over 
the whole IM-MI band by means of the system param¬ 
eters, which brings to the idea of wideband parametric 
amplification in multimode fibers. 


IV. A PHYSICAL INSIGHT INTO THE IM-MI 

The physical model resumed by Eq.([8]) offers a semi- 
analytical solution that can precisely quantify the IM-MI 
growth in terms of eigenvalues and eigenvectors of the 
matrix M. 

A natural question that arises is how to control the IM- 
MI gain bands by means of the system parameters. From 
this point of view the model resumed by Eq.([8]), although 
accurate, lacks simplicity as it makes necessary to com¬ 
pute at each detuning f2 the eigenvalues of M. 

For this reason in this section we aim at gaining a deeper 
physical insight by looking for some basic model that can 
easily explain the main features of the IM-MI process. 
As outlined in the Introduction, a simple way of thinking 
the IM-MI in a multimode fiber is to consider separately 
the nonlinear interactions between all the possible cou- 


















0 I —-- . rL — v- a v - 

-1.1 -1 -0.9 


Detuning v 


FIG. 4. Amplitude amplification factor for the 2x— mode 
when L = 5. Frequencies around v = 0 (where the pump is lo¬ 
cated) have been filtered out. Panel (a): band from v = 0.85 
to v = 1.1; (b): band from v = 0 to v = 0.65; (c): band from 
v = —0.65 to v = 0; (d) band from v = —1.1 to v = —0.85. 
In gray: A^, x obtained by numerical solution of Eq.|3j) by 
split-step method. In red: the estimation A 2 X ,est obtained by 
Eq.(13l. The vertical dashed lines in panel (c) delimit the 
band of Bg centered at u = —0.43. The inset in panel (c) 
shows a zoom of A 2 X around v = —0.43; the gain curve Bq 
(blue) is also reported. The difference between the peaks of 
the two curves is nearly 0.07, that is to say A 2 X ~ Bg — 0.07 
in the band of Bg- 


pies of spatial and polarization modes. This approach, 
which we call bimodal-MI model, has been adopted in 
previous works concerning the IM-MI mm- In this 
section and in the next three sections we revise this basic 
model in order to find some useful analytical formulas for 
the IM-MI gain as well as to discuss the impact of the 
TOD and losses in the IM-MI dynamics. 

If one takes into account for the interaction between the 
na- mode and the mb- mode, then Eq.Q should be rewrit¬ 



FIG. 5. Same as in Fig. |4](c), but when L = 16. The inset 
shows a zoom of A 2 X around v = —0.43; the gain curve Bf 
(magenta) is also reported. The difference between the two 
functions is nearly 0.21, that is to say A 2 X ~ Bf — 0.21 in the 
band of Bg- 



FIG. 6. Amplitude amplification factor for the 4x —mode 
when L = 5 and in the band from v = —0.65 to v = 
0 .Frequencies around v = 0 (where the pump is located) have 
been filtered out. In gray: A± x obtained by numerical solu¬ 
tion of Eq. <§ by split-step method. In red: the estimation 
A- 4 X ,est obtained by Eq.(131. 


ten by neglecting all the terms except for those involving 
the sidebands 


ds ri 


dz 


— ^(0na(~) 4” bsCnn\Pna\ )'5na - t _ 


ib s C nn p 2 na i* na exp(-iApMz - iA/3^ x) z)+ 

ibmnCmnPmbPnaSmbexp{iA^ S) Z - iAp^z) + 

'ib mn C rnn p mb p na i* mb exp(-iA/3 < ^j } l) z - iA^ s) z) 


(15) 


where b mn = 6|| if a = b and n/m; b mn = b± if a ^ b 
and n/ra; b mn = bx if a ^ b and n = m; b mn = 0 if 
a = b and n = m. An equation similar to Eq. (15) could 
be written for i na by exchanging s F4 i , for s mb by ex¬ 
changing n to and a -f4 6, and for i mb by exchanging 
s F4 i, n F4 m and a F4 b. 


The first term in the right-hand-side of Eq. (15) is here 


called phase-term as it is responsible for the phase mod¬ 
ulation of s na . The second term is called self-MI-term as 
it accounts for the degenerate four-wave mixing (FWM) 
process where two na— pump photons are converted to 
an na—idler photon and a na— signal photon, leading 












































9 


to their amplification. Similarly, the last term is called 
cross-MI-term as it accounts for the FWM where a 
na —pump photon and a mb —pump photon are converted 
to an na —signal photon and a mb— idler photon. The 
third term is called not-MI-term as it does not account 
for sideband amplification. 

We point out that if a = b and n = m we finally get a 
system of 2 equations for s na and i na that closely recall 
the system of equations describing the scalar MI (SMI) 
in a single-mode fiber [ 21 ] ; here we thus refer to this in¬ 
stance as SMI na . 

Otherwise, if a 0 b and/or n / m we get a system of 4 
equations for s na , s m b , i n .a and i m b that recall the sys¬ 
tem of equations describing the vectorial MI (VMI) in 
single-mode fibers m ■ Here we refer to this instance as 
MI na -mb- 

By generalizing the SMI in single-mode fibers [21] to the 
most general case of the SMI na in a multimode fiber, we 
may easily infer that the gain curve B na related to the 
SMI na process reads as : 


B na = -l/Val^c -^ 2 )!/ 2 


v \P2,nn \ ' 


1/2 


(16) 


which is valid under the anomalous dispersion condi¬ 
tion /? 2 na < 0 and where fl c is the upper cut-off frequency. 
As regards the MI na - m b process, a rich and complex dy¬ 
namics is found that strictly depends on the ratio bs/b mn 
(see for example [IB! for the VMI in single-mode Hi-Bi 
fibers). The general way of handling the system of 4 equa¬ 
tions for s na , s m b , ina and i m b consists in transforming it 
in an eigenvalue equation by means of a transformation 
similar to Eq. 0- and then to compute the corresponding 
dispersion relation. The problem is that, unless peculiar 
instances, the dispersion relation is a fourth-order poly¬ 
nomial equation for Q requiring numerical computation 
to be solved. 

On the other hand, different spatial modes typically ex¬ 
hibit a large group velocity mismatch (GVM) and the 
same occurs to the polarization modes in a Hi-Bi fiber. 
In this case the phase-matching conditions for the self- 
MI term, for the cross-MI-term and for the not-MI-term 


in Eq. (15) are reached at largely different detunings. For 


this reason, with the aim of studying the amplification 
induced by the cross-MI-term, we keep it along with the 
phase-term in Eq. (15), whereas we neglect the self-MI- 


term and the not-MI-term. Doing so, we can derive the 
following analytical estimate for the gain B na _ m b related 
to the MIna—mb process (see Appendix II for details): 


Bna—mb — (&i — — Slpp) 2 /^PK,Nl) 

flpK = ^PK,L + ^PK,NL 

n pk,L = |(2 Dps )- 1 (-02 ± {02 2 - 41^3 A,) 1/2 ) I 
^PK,NL = 2A;2/(3 ^l 2 p KL Bp 3 + 2Qpk,l02 + D v ) 

kl “ b m nC'mn\Pna\\Pmb\ 

h = (b S /2){C nn \Pna\ 2 + C mm \p mb \ 2 ) (17) 


where Upp stands for the peak-gain detuning that 
here is written as the sum between a linear contribution 
flpK,L and a nonlinear contribution f Ipk,nl] 02 is 
the average GVD {02,na + 02,mb)/2; Dp 3 indicates 
the TOD difference { 0 3na — /? 3 mb )/6 and D v is the 
GVM = v^-v~l ’ ’ 

The ± operator in the third line of Eq. ( 1 T 7 J) reveals 
that two distinct gain bands are generally found when 
Dp 3 ^ 0, let us say This issue is carefully 

addressed in the next section and clearly shows that the 
TOD, which does not affect MI in single-mode fibers, 
can play an important role in the IM-MI dynamics 
instead. 

In Fig. [ 7 ] the gain curves B na and B na _ m b, computed 
according to the analytical estimations Eqs.( [l 6 j [Tt| , 
are compared to the gain curves Ba...j of the full model 
discussed in the previous section and displayed in Fig. [2] 
The two gain curves b[ 1 J_ 2x and B^_ 2x clearly match 
with Be and Bh , respectively. Similarly b[ x _ 4x and 
b[ 2 J_ 4x match respectively with Be and Ba ; B^J_ 3x 

and b!0J_ 3x match respectively with Bq and Bp, B 3 J_ 4x 
( 2 ) 

and B 3 J_ 4x match respectively with the right-side of Bp 
and with Bp. 


The SMI 2x and the SMI 4x processes occur because 
of the anomalous dispersion conditions 02 , 2 x < 0 an d 
02,Ax < 0 (see Table I). Moreover, their corresponding 
bands B 2 X and B^ x match with Bj and the left-side of 
Bp, respectively. 

We note the absence of the two gains related to the 

Ml 2 x- 4 x process, which is due to the large TOD dif- 
- 2 

ference that makes 02 — 4Dp 3 D v < 0 and thus f 1pk,l 
in Eq.( |17| ) not real-valued. As explained in next section, 
these bands may appear by introducing higher-order dis¬ 
persion terms that we have not contemplated for the sake 
of simplicity. 

The first of the two gain curves related to the MI\ x _ 3x 
process, that is , matches with Bp- The second 

is centered at v = 4.6, out of the frequency window from 
v — —1.1 to v = 1.1 displayed in Fig. [2] and thus it is 
not represented. 

It can be seen that the agreement between the approx¬ 
imated curves B na , B na -mb and the gain curves Ba...j 
is quite good at all detunings. This outcome confirms 
that we can primarily ascribe each one of the IM-MI gain 
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FIG. 7. Comparison between the gain curves Ba...j computed 
from the full model Eq.|8| (blue solid lines) and the analytical 
estimations calculated from Eq.(|16|[l7|) (black and red circles). 


bands to one of the possible SMI na and MI na - m b pro¬ 
cesses in the fiber. As a consequence, we can employ 
Eqs.( 16 17) to get an analytical estimation of the IM- 
MI gain. 

Yet, we stress that some discrepancy is present between 
the approximated gain curves and Ba...j, which is due 
to having neglected the self-MI-term and/or the not-MI 
term in Eq. (15) as well as the interplay between all the 
fiber modes. 

Indeed Eq.( 15) describes the dynamics of 1 (if n = m and 
a = b) or 2 (if n ^ m or a ^ b ) modes that we call phase- 
matched (PM) modes and that are amplified thanks to a 
phase-matched FWM , which is the self-MI-term or the 
cross-MI-term of Eq.(fT5|). On the other hand, Eq.(15) 


totally disregards the other modes, that we call not-PM 
modes. 

As already stressed, the not-PM modes also undergo am¬ 
plification and their dynamics can be described only by 
means of the full model Eq.®. Roughly speaking, the 


not-PM modes are amplified because they are nonlinearly 
coupled to the PM modes by means of a not-PM FWM. 
Let us consider for example the gain band Be, which 


is related to the MI^_ ix process and where the side¬ 
bands of the lx—mode and the 4x—mode are amplified 
by means of a PM FWM process. The 2x—mode does 
note take part to this PM FWM (it is a not-PM mode 
in the band of Be) but even so it is nonlinearly cou¬ 
pled to the lx—mode and the 4x—mode: indeed from 
Eq.([6]), written for the 2x—signal, we find that d z s 2x 


depends on p^p^i^expt—iAp^)^ z — iA[3^ s> z) and 

PixP2xi* x exp{-iA/3^ z - iA^^ z) . Although these 
last two FWM terms are not phase-matched in the band 
of Be, still they lead to the amplification of s 2x because 
the idler i\ x and ii x grow exponentially. 

We can therefore conclude that, although the basic model 
of Eq. (15) sheds an important light on the IM-MI dy¬ 
namics and allows deducing some useful analytical ap¬ 
proximation of the IM-MI gain curves, however it cannot 
describe the MI growth of the not-PM modes. At this 
purpose, the full-model Eq.([8]) should be exploited, as it 
provides both the eigenvalues and the eigenvectors in or¬ 
der to correctly characterize the spatial evolution of each 
mode. Furthermore, it is only by means of the full-model 
that we can precisely calculate the IM-MI gain curves 
without any approximation. 


□ (p.s). 


V. IMPACT OF THIRD-ORDER DISPERSION: 
SECONDARY MI BANDS 


Typically the large GVM between the na —mode and 
the mb —mode in the M na _ m b process shifts the corre¬ 
sponding gain B na _ m b towards high frequencies; conse¬ 
quently, the role of higher-order dispersion terms should 
be carefully analyzed. At the best of our knowledge, 
this has not been done in previous works concerning the 
bimodal-MI jTTUTB] . However, a full treatment is com¬ 
plex and out of the scope of this paper, therefore in this 
section we just provide an outline of the influence of TOD 
on the IM-MI dynamics. 

Equation reveals that two distinct bands are gener¬ 
ally found when \Dp 3 \ > 0 and it thus gives evidence of 
the remarkable impact of TOD on the IM-MI. 

A similar behavior has been observed by Nithyanandan 
et al. in [[22] in the case of two pumps co-propagating 
at different frequencies in a single-mode fiber. In that 
case a large TOD difference leading to the formation of 
2 distinct bands was induced by a large frequency detun¬ 
ing between the two pumps, whereas in our case (single¬ 
pump) it originates from the different dispersion charac¬ 
teristics of different modes. 

In order to better understand the role played by TOD, 
it is instructive to explore the limit \Dp 3 \ -A 0. From 
Eq.p7| we see that in this limit the 2 bands are roughly 
centered around the peaks ss \D v /f3 2 \ and Opj^- ~ 
\/3 2 /1, respectively (we have neglected the nonlinear 
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contribution £Ipk,nl) ■ 

The first gain band, centered around is practically 
unaffected by the TOD. It is the conventional band dis¬ 
cussed in previous works concerning the bimodal-MI. 

( 2 ) 

The second band,centered around Q PK , is strongly de¬ 
pendent on Dp 3 and represents a secondary MI spectrum 
in addition to the conventional one. It is typically located 
far away from the pump and for this reason higher-order 
dispersion terms besides the TOD could be not negligi¬ 
ble. 

According to Eq. the more \Dp 3 | increases the more 
the two bands move closer each other, but in a different 
way depending on the sign of Dp 3 and D v (see Fig. [8|. 
If sign(Dp 3 ) = —sign(D v ) then both the bands move to¬ 
wards lower frequencies. Otherwise the conventional one 
moves towards higher frequencies and the secondary one 
towards lower frequencies, until they theoretically over¬ 
lap once that /3| = ADp 3 D v . On the other hand, in the 
limit /3 2 —> 4:DpsD v the higher order terms besides the 
TOD play a not-negligible role, so that their influence 
on the conventional and secondary band should be taken 
into account. 
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FIG. 8. Schematic representation of TOD impact. Panel (a): 
conventional and secondary bands related to the MI na - m b 
process in the limit | Dp 3 \ —>• 0. Panel(b): case of sign(Dp 3 ) = 
sign(D v )\ for increasing values of | Dp 3 | the conventional band 
blue-shifts, whereas the secondary band red-shifts. Panel(c): 
case of sign(Dp 3 ) = —sign(D v ): for increasing values of \Dp 3 \ 
both the bands red-shift. 


VI. IMPACT OF THE PUMP MODAL POWER 
DISTRIBUTION 

As the higher-order dispersion terms, so the modal 
power distribution of the pump plays an important role 
in the IM-MI dynamics. This is evident considering that 
different power distributions produce different eigenvec¬ 
tors and eigenvalues of M. 

In this section we want to focus on the impact of the 
pump modal powers on the IM-MI gain bands. 


Past works concerning the bimodal-MI mm have not 
put in evidence the nonlinear contribution f ^pk,nl to 
the peak-gain detuning flpK related to the MI na _ m j J 
process. Although generally Qpk,nl « U pk , however 
^pk,nl has the same order of magnitude of the band¬ 
width of Bna-mjj and should therefore taken into account 
if we want to precisely locate this gain band. At this pur¬ 
pose we point out that, according to Eq. 0 , B na — rn i ) as 
function of D is a semi-ellipse centered in f2 pk and whose 
bandwidth is 2 flpK,NL(ki/k 2 )- 

It is also interesting to note that the peak-gain k\ = 
bmnCmn\Pna\\Pmb\ of -B na -m b is directly controlled by 
the two pump modal powers |p no | 2 and \p m b\ 2 ■ Under 
the constraint |pna| 2 + \Pmb\ 2 = constant , k\ is max¬ 
imized when |p no | 2 = \pmb 1 2 , namely, when the pump 
power is equally distributed over the na —mode and the 
mb— mode. 

Contrary to the MI na _ m b, the SMI na process occurs 
without any GVM or TOD contribution to the phase¬ 
matching condition of the self-Mi-term in Eq. (15). In¬ 
deed, phase-matching is achieved through direct compen¬ 
sation of second-order dispersion by self-focusing nonlin¬ 
earity and as consequence both the position of the peak- 
gain detuning and the value of the peak-gain are strongly 
dependent on the pump power \p na \ 2 - 


VII. IMPACT OF MODAL PROPAGATION 
LOSSES 


Typically higher-order spatial modes suffer larger at¬ 
tenuation than lower order modes, which makes inter¬ 
esting to study the influence of losses in the context of 
IM-MI. 

Past works concerning the MI in single mode lossy fibers 
have demonstrated that the MI dynamics could be signifi¬ 
cantly influenced by fiber losses f [23U25] ): indeed an opti¬ 
mum propagation length exists for which the MI growth 
is maximized and beyond which it is gradually annihi¬ 
lated by losses, and the MI spectrum experiences a con¬ 
tinual frequency shifts towards the pump. 

In this section for simplicity of notation we assume that 
losses are not dependent on SI; however S7—dependent 
losses could be easily treated. 

Modal losses are included by adding the terms ia na 
and —ia na to the elements M sajSO [n, nland M JQj i Q [n, n] 
of the matrix M defined in Section ITTl respectively ( 
1 < n < N,a = {x,y} ). In addition, losses cause 
an exponential decay of the pump amplitudes, therefore 
the terms |p na | of the matrix M should be replaced by 
| p na (z = 0)\exp(—a na z). In this way the coefficients of 
the matrix M are no longer constant, which prevents the 
existence of a simple solution to the eigenvalue problem 
d z v = iMv of Eq.([8]). 

On the other hand we are interested to the regime in 
which the z-dependent coefficients \p na {0)\exp(—a na z) 
are slowly decaying if compared to the fast Mi-growth 
of the sidebands, otherwise a fiber-amplifier would be 
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low-efficient. In this regime the solution v could be well 
approximated by a Magnus series expansion truncated 
at the first order ([2B]), that is v(L) = ea;p(iM av L)v(0), 
with M av = L~ x Jq M (z)dz. 

This means that when losses are taken into account then 
the vector v(L) is still obtained by Eq.(10) provided that 
we compute the eigenvalues Xk and the eigenvectors Wk 
of the matrix M av , which is constructed from the ma¬ 
trix M defined in Section [IT] by replacing the terms \p n a\ 
with their spatial average L~ x \p na (0)\exp(—a na z)dz 
and by adding the loss coefficients ia na and —ia na as 
previously explained. 

Note that large losses may prevent the existence 
of MI phenomena: in a single-mode fiber, this 
corresponds to the case in which propagation 
losses are larger than the MI gain. We should 
therefore wonder if a given input pump power dis¬ 
tribution allows or not for IM-MI. The response 
comes from the eigenvalues of the modified ma¬ 
trix M av evaluated at the fiber entry (i.e. fixing 
L = 0): if at least one eigenvalue with negative 
imaginary part exist, then IM-MI will occur. 

In order to assess the validity of our approach, we repeat 
the simulation described in Section |III| but including the 
following amplitude losses (per nonlinear length L nl ^)\ 
a lx = 0.005; a 2x = 0.01 ; a 3x = 0.015; a 4x = 0.02. 
In the case under analysis such losses are much larger 
than typical losses in silica fibers at telecommunication 
wavelengths^^, = 0.005 corresponds to a real power 
loss of about 0.43 dB/m),nonetheless they are small com¬ 
pared to the IM-MI gain; moreover, here our purpose is 
to show the robustness of our approach even in presence 
of a relatively fast pump decaying. 

In Fig. [9] we show the good agreement between the am¬ 
plitude amplification factor A 2x , obtained by nume rica l 
simulations, and its estimate A 2xest , obtained by Eq.( 13) 
where the dominant gain and the corresponding eigenvec¬ 
tor are those associated to the matrix M av . 

In this example losses have a deep impact over the MI 
growth of the 2x— mode. We have seen in Section |III| that 
in absence of losses the gain g 2 = 0.9 is larger than g\ = 
0.71, which makes g 2 to be dominant, for the 2x— mode, 
starting from L = 15.8. For this reason, in L = 16 
and around v = —0.43, A 2x ss Bp + L~ x Log(\w 2 [2]\) 
(see dashed blue line in the band of Bq in Fig. [ 9 ]). On 
the other hand when losses are included, we find that 
for L = 16 the new gains g 2 and g 1 , calculated from 
M av , are respectively g 2 = 0.45 and g 4 = 0.56, which 
makes gi to be dominant instead of g 2 and therefore 
A 2x ~ B g + L~ 1 Log( |iDi[2]|) (Fig. [ 9 ) gray line). 
Following a treatment similar to that exposed in Ap¬ 
pendix II the estimate of B na _ m b given by Eq.([T7|) could 
be recalculated by taking into account for losses. The 
matrix M' of Appendix II should be rewritten by adding 
the loss coefficients and by replacing the pump ampli¬ 
tudes with their corresponding spatial averages. 

After some algebra we find that the peak-gain 


Bjia—mb pk ma y be approximated as: 


B 


(Is) 

na—mb,PK 


= B na - mb)PK (l - exp(-2aL) s j / (2aL) - , 


(18) 


where the superscript (Is) stands for losses, whereas 
B na -mb,PK is the peak-gain in absence of losses and a in¬ 
dicates the average ( a na +a m b)/2. Interestingly enough, 
the drop of the peak-gain is thus related to the average 
between the loss coefficient of the na —mode and the loss 
coefficient of the mb— mode. 

A similar equation holds true for the peak-gain B y na PK of 
the the SMI na process, after replacement of B na _ m b,pK 
with B na px and of a with a na in Eq. 18 


Furthermore, the peak-gain detuning related to the 
SMI na undergoes a red-shift, which is completely simi¬ 
lar to the MI frequency drift already predicted in single¬ 
mode fibers and which can be estimated as follows: 


= ttpK (l ^ exp(-a na L )) / (a na L) (19) 

On the contrary, in the processes only the 

nonlinear contribution pk,nl is shifted by losses, but 
not the linear contribution &pk,l which is significantly 
larger than f ^pk,nl- As a consequence,the global peak 
shift induced by losses is typically negligible. 



FIG. 9. Amplitude amplification factor for the 2a:—mode 
when L = 16 and in the band from v = —0.65 to v = 0. In 
gray: A 2x obtained by numerical solution of Eq. ([ 3 J) by split- 
step method when propagation losses are included (aia, = 
0.005; a 2x = 0.01 ; 0 : 3 a, = 0.015; a 4x = 0.02). In red: the es¬ 
timation A 2Xj est obtained by Eq. (13) using eigenvectors and 
eigenvalues of M av . In dashed blue: A 2x obtained by numer¬ 
ical solution of Eq. CD but with zero losses. 


VIII. LIMITS OF VALIDITY 

In this section we discuss the limits of validity of the 
main results deduced in this paper. 

We remind that in Eq.([2]) several nonlinear terms 
have been neglected because considered fast-oscillating. 
As an example, let us take the nonlinear term 

CnnmndPrixB r nx^'l,nnmn: with Tl ^ TYl and Cnnmn 7 ^ 
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0, whose corresponding wavevector mismatch is 
A/?i iTmmn = /3 nx (u> p ) - /3 mx (uj p ). Typically the propa¬ 
gation constants /3 nx (ui) and f3 mx (ui) are largely different 
whatever the frequency u> is, which justifies the statement 
\^‘0\,nnmnL nl\ » 0. On the other hand, for a proper 
fiber design these 2 propagation constants could become 
equal at a peculiar frequency, let us say u. If this is the 
case and the pump frequency w p coincides with u>, then 
this nonlinear term must be taken into account in Eq.|2]). 
What said is generally true for anyone of the terms that 
we have neglected. 

Note also that, although the neglected terms are typi¬ 
cally truly fast oscillating, however they could lead to 
some MI process that we have not accounted in our 
analysis. As an example, let us consider once again 
the term c nnmn Pf x P* lx e\ jnnmn previously discussed, as¬ 
suming / w so that it is effectively fast-oscillating. 
Its decomposition by means of Eq. @ would lead to the 
presence of the term c nn mnPn X i* mx exp{i2l3 p ^ nx - i/3 itrnx - 
iPs,nx) in Eq. 0, which is associated to a FWM process 
where energy is transferred from the nx —pump to the 
nx —signal and the mx— idler. Once again, these consid¬ 
erations are generally true for anyone of the terms that 
we have dismissed in Eq.([2]); consequently, a group of MI 
bands exist that we have ignored in our analysis. How¬ 
ever, these MI bands are typically located far away from 
the gain bands predicted by our theory; moreover the 
corresponding FWM processes are very sensitive to small 
fiber imperfections so that their typical coherence length 
is less than 1 meter (ini), which practically prevent an 
efficient MI growth in fibers longer than a few meters. 
For this reason we can safely assume that these neglected 
Mi-bands do not affect the analysis developed in this 
paper. On the other hand, for very high-input powers 
and/or a large number of propagating modes in short 
fibers, some of the IM-MI gain bands predicted by our 
analysis would be located far away from the pump and 
could thus overlap with the neglected Mi-bands, which 
may invalidate our results in the portion of the spectrum 
where overlap occurs. 

Another important issue that we have not considered in 
our analysis is the effect of polarization mode dispersion 
(PMD) as well as of spatial mode dispersion (SMD) re¬ 
lated to groups of quasi-degenerate modes, which could 
greatly impair the IM-MI growth in telecom fibers m- 
Indeed, in the case of telecom fibers, Eq.Q is strictly 
valid in the regime of zero-PMD and zero-SMD only. 

A proper analysis taking into account for the PMD and 
SMD could be done following a treatment similar to 
that proposed in m in the case of single-mode fibers. 
The final result would be the presence of exponentially- 
decaying nonlinear coefficients in Eq.([6]), which is com¬ 
pletely analogue to the presence of exponentially decay¬ 
ing pump amplitudes when propagation losses are taken 
into account. We could therefore employ the Magnus ex¬ 
pansion proposed in Section |VII| to deal with both losses 
and PMD/SMD. On the other hand, the derivation of the 
exponentially-decaying nonlinear coefficients in presence 


of PMD and SMD is complex and out of the scope of this 
paper. 


IX. CONCLUSIONS AND PERSPECTIVES 


In this paper we have presented a detailed theoretical 
and numerical analysis of the IM-MI in multimode 
fibers, which brings to the idea of wideband multimode 
parametric amplification. We can summarize the main 
results as follows: 

-Eigenvalue equation for multimode fibers- One of 
the main outcomes of this work consists in the lineariza¬ 
tion of the coupled nonlinear Schrodinger equations 
describing the propagation in a multimode fiber and in 
their following transformation in an eigenvalue problem. 
In this way the modal amplification can be described by 
means of the eigenvectors and eigenvalues of a matrix 
M that contains the information about the dispersion 
characteristics of the modes and the modal power dis¬ 
tribution of the pump. This result generalizes previous 
studies of MI with 2 spatial or polarization modes to 
the most generate case of N > 2 interacting modes. The 
computation of the eigenvectors and eigenvalues of M 
as function of the pump-sideband detuning S2 permits 
to completely characterize the IM-MI. In particular, the 
eigenvalues with negative imaginary part, called as usual 
gains, lead to sideband amplification. 

-Dominant gain and wideband amplification- For a 
given detuning several gains could be found, each one 
playing the role of dominant gain at different fiber posi¬ 
tions and for different fiber modes. The amplification of 
each mode is thus controlled by its dominant gain as well 
as the corresponding eigenvector, which in turn depend 
on the system parameters. This issue is well highlighted 
by Eq. which offers an analytical estimate of the 
modal amplification. Most importantly, as highlighted 
in Section [TlJ each mode undergoes amplification in the 
whole IM-MI band. 


-Physical interpretation of the multimode amplification- 
Although the eigenvalues and eigenvectors of M can 
completely characterize the multimode amplification, 
however they do not provide a simple and intuitive 
picture of the physics behind the IM-MI. At this purpose 
we considered separately the nonlinear interactions 
between all the possible couples of modes, which allows 
decomposing the complex IM-MI dynamics in a set of 
bimodal-MI processes. The bimodal nonlinear interac¬ 
tion between the na —mode and the mb —mode gives rise 
to a MI process, called MI na _ mb , completely analogous 
to the vectorial MI in single-mode fibers. Similarly the 
self nonlinear interaction of the na —mode gives rise to 
the SMI na process, which is analogous to the scalar MI 
in single-mode fibers. These analogies permit to find 
out the analytical estimates Eqs.( 16 17) approximating 
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the IM-MI gain and to interpret the wide IM-MI band 
as the union of the bands related to each MI na - m b and 
SMI na process. As a result, the useful bandwidth is in 
principle much larger than in the case of single-mode 
fibers. 

-Phase-matched and not-phase-matched modes- In 
the simple bimodal model describing the MI na - m b and 
the SMI na processes the na —mode and the mb— mode 
play the role of phase-matched modes, namely, they 
are amplified by means of a phase-matched FWM. 
Nevertheless, this simple model cannot give any infor¬ 
mation about the other modes, which we call not-phase 
matched and which undergo amplification because they 
are nonlinearly coupled to the phase-matched modes. 
Therefore, although it sheds an important light on the 
IM-MI dynamics, the bimodal model does not offer a 
complete characterization of the IM-MI amplification, 
which is provided by the eigenvectors and the eigenvalues 
of M instead. 

-Influence of TOD and losses- We have concluded 
our work by analyzing the impact of TOD and of losses 
on the IM-MI dynamics. We have put in evidence that 
the large GVM between different spatial modes may shift 
the IM-MI gain bands towards high frequencies, so that 
higher-order dispersion terms should be contemplated 
in order to correctly describe the IM-MI. Differently 
from single-mode fibers, the TOD difference between 
different modes plays an important role and leads to 
the formation of secondary Mi-bands in addition to the 
conventional ones. 

Contrary to the GVM, propagation losses may shift the 
gain bands towards low frequencies. Their impact (gain 
drop and frequency shift) can be evaluated by properly 
modifying the matrix M and resorting to a Magnus 
expansion which accounts for the spatial decaying of the 
pump amplitudes. 

-Future perspectives- The outcomes exposed in this 
paper pave the way towards the implementation of 
wideband multimode parametric amplifiers where the 
amplification of each mode could be selectively con¬ 
trolled by means of the system parameters. For this 
reason a natural evolution of this work is the study 
of the optimization of the pump modal distribution 
and of the fiber parameters in order to maximize the 
amplification of one or more modes in a desired band. A 
further important issue to address will concern 
the system scalability. In addition, although limited 
to the case of single-core step-index fibers, the theory 
here developed could be conveniently modified so as to 
describe IM-MI phenomena in multicore fibers or pho¬ 
tonic crystal fibers, which provide a unique opportunity 
for tailoring the modal dispersion characteristics and 
the nonlinear coupling coefficients. Finally, this work 
may find a useful application in the description of 
the complex spatio-temportal soliton dynamics 


in multimode fibers, which recently has been 
widely discussed mmm and represents a 
truly hot-topic in optical-physics. 
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APPENDIX 1 


In this Appendix we derive Eq. (13), which represents 


an analytical estimate for the amplification factor 
related to a weak input background noise amplified by 
the IM-MI process. 

Each input modal field is the sum between a pump 
component p no (0,f) and a noise r na (0,t). For the sake 
of clarity in this Appendix we refer to the lx—mode but 
what follows can be easily generalized to anyone of the 
propagating modes. 

The power spectrum Ri x (z,u) = \F{r\ x {z, f)}| 2 {F 
indicates the Fourier transform) typically exhibits strong 
and fast fluctuations due to the random nature of r\ x 


(Fig. Al-1). In order to correctly estimate the 
noise spectrum, and thus to reduce the spectral 
fluctuations, a statistical average should be per¬ 
formed over several realizations of input noise. 
Clearly, this approach is extremely time consum¬ 
ing. An alternative option, which turns out to be 
equivalent to the aforementioned approach, con¬ 
sists in averaging the power spectrum in a narrow 
band of width b where several fluctuations are 
included, so that the averaged power spectr um re ads 

as Rl x ,avg(z, v) = 6 _1 J ^/2 R\x(z, V’W (Fig. 


Al-11. 



0 0.2 0.4 0.6 0.8 1 

Detuning v 


FIG. Al-1. Comparison between R\ x (red) and Ri Xl av g 
(black). 

This averaging process is practically equivalent 
to applying a moving-average filter to Ri x (z, v)\ it 
allows suppressing the fast fluctuations in the spec¬ 
trum and recovering the real noise power level around a 
certain frequency. 

In order to find an estimate of the amplification factor 
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A± x , we start from Eq.(10) written for the la;—mode, 
namely: 


AN 

v[l }{z,is) = ^ ~2ckWk[l]exp(i\kZ ) (Al-1) 

k =1 


where v[l]{z,is) indicates the amplitude of the 
lx—mode at the position 2 and at the detuning fre¬ 
quency v, that is |u[l] ( 2 , v) | 2 = R\ x {z,v), and the 
wavevectors are normalized so that |wfc • w^l = 1. 
For notational simplicity, in this Appendix we omit to 
indicate the dependence of Cfc, Wk and A k on i/, although 
they are actually ^—dependent; indeed Wk(V) and A kiy) 
(1 < fc < AN) are respectively the eigenvectors and the 
eigenvalues of the matrix M computed at the frequency 
v, whereas Cfe(^) = (v(0, v) • Dwk)/(wk • Dwk). We 
also omit to indicate the interval of integration, so we 
indicate simply with f. 

The elements x[n](0, u) and v[n + ./V](0, v) of v(0, v) 
(1 < n < N) represent the input noise amplitude 
for the nx— mode and the ny— mode at the Stokes 
frequency u, respectively, whereas v[n + 2iV](0, v) 
and v[n + 3iV] (0, z^) are the noise amplitudes for 
the nx—mode and ny—mode at the anti-Stokes 
frequency — u\ therefore |u[n](0, v )\ 2 = R n x{ 0, i'), 

|x[n + AT](0, v)\ 2 = R ny { 0,j/), | v[n + 2iV](0, v)\ 2 = 
Rnx(0,~v), \v[n + 3iV](0, p)\ 2 = R ny {0, —v). Further¬ 
more, the integrals b~ l J |v[n](0, =R nx ,avg{ 0, z') 

and b~ 1 f \v[n + N](0,ijj)\ 2 dilj =R n y,avg{ 0, i/) indicate re¬ 
spectively the input noise power related to the nx—mode 
and to the ny—mode in the Stokes band [u — 6/2, v + b/ 2]; 
similarly 5 _1 f \v[n + 2V](0,■0)| 2 ch/ ) =R n x,avg( 0,~v) 
and 5 _1 J'\v[n + 3iV](0, i/j)\ 2 dip =R ny , aV g(0,-v) indi¬ 
cate the input noise powers in the anti-Stokes band 
[-v - 5/2, -v + 5/2]. 

Since the input background noises are white and inde¬ 
pendent of each other, we can safely assume that the 
averaged spectra R nx ,av g { 0, n) and R ny ,av g {0, v) are 
independent of the polarization, of n and v, so that 
Rna^avg (0, Av)ttR\ x ^ aV g (0, v) (Cl — {x,y},l ^ Tl A N ) 
and therefore / |x[fc](0, i/)| 2 ch/>~/ |x[l](0, ij))\ 2 d^ = 


bRix,avg{0, v) {Wk : 1 < k < AN). 

We may approximate Eq.(Al-l) by taking into ac¬ 
count only the dominant gain gdom.i for x[l] at the 
position 2 and the associated eigenvector w dom,ii 
that is v[l)(z,v) « Cd om ,iw d om,i [l}exp{g do m,iz). 

In this way we could rewrite Rix,av g {z,v) = 
b^ 1 f \v[l}(z,i/j)\ 2 di/j=b- 1 p 1 f \c do m,i\ 2 d0>, where 

Pi = \wdom,i[]]\ 2 exp{2gdom,iz) and we have assumed 
Wdom ,i [1] and gdom ,l to be constant in the small band 
[v — 5/2, v + 5/2] so that they can be taken outside the 
integral. 

When calculating f \cdom,i\ 2 di/j we should take 
into account that couples of different input modal 
noises are independent each other, therefore 
f x[n](0, v) ■ v[m]* (0, v)diji k 0 if n / m and then 


/ \cdom,l\ 2 d-lll « P2 Y?k=l \ W dom, 1 [A] | 2 / \v[k\(0, 1p)\ 2 d'lfj , 

with p 2 = |w dom ,i • Dw domj i| -2 . 

Being J2i=i\ w dom,i[k}\ 2 = \™dom,i • w dor71jl | = 1 

and f \v[k](0,ip)\ 2 d'i/j « bR 1X:avg (0, v), we con¬ 

clude that f \cdom,i\ 2 dv=p 2 bR 1Xi av g (0,v), and then 
Rix,avg{z, v) = PiP2Rix,avg(0, v), from which Eq. (J13J) is 
easily derived. 


APPENDIX 2 


In this Appendix we derive Eq.( |17| ), which represents 
an estimate for the gain B na _ m b related to the MI na _ m b 
process. 


As suggested in Section IV in Eq.(|l5[) we may just re¬ 


tain the phase-term and the cross-MI-term. As a re¬ 
sult, the system of 4 CNLSE describing the dynamics of 
SnaAna, s m b,imb is splitted in two distinct subsystems: 
the first accounting for the interaction between s na and 
i m b ; the second for the interaction between s m b and i na . 
We indicate with MI sna _i m b and MIi na _ sm b the corre¬ 
sponding MI processes. 

As regards the MI sna _ im b : making use of Eqs.( [5j [7]) 
we finally get the eigenvalue problem d z V = iM'v , be¬ 
ing v' = [s na imb\ T an( f 2x2 matrix whose ele¬ 

ments are respectively M'[l, 1] = M sajSa [n, n], M'[2, 2] = 
-M ib ,ib [m,m] and M'[l,2] = -M'[2,1] = M sa , ib [n,m] 
(matrix M is defined in Section 0- 
The gain B na _ m b is given by the eigenvalues A of M' 
with negative imaginary part, namely: 


zb = (M'[l, 2] 2 - 4 _1 (M'[1,1] - M'[2, 2]) 2 ) 


1/2 


(A2-1) 


According to Eq.(A2-l) the peak-gain M'[l,2] = 
bmnCmn\Pna\\Pmb\ is reached at the peak-gain detuning 
f Ipk for which the phase-matching condition M'[l, 1] — 
M'[2,2] = 0 occurs. The term M'[l,l] — M'[2,2] = 
L + NL is the sum of a linear part L and a nonlinear 
part NL that read as: 


NL — bs{C nn \Pna |~ 4” Cjjim |Pm£>|~) 

L(n) = + A/3^ = Dp 3 n 3 + p 2 n 2 + D v n 

(A2-2) 


where A/?io’ s) and A/3^ are given by Eq.(lT I, D t 33 = 
(/^3,na ^3 ,mb)/b, 02 — {02 ,na 4” 02,mb )and D v = 

v na — v mb- Note that according to this notation Ll is 
positive-valued. In this Appendix we assume that the 
coupling coefficients C nn and C mrn are independent of Ll 
in the band of B na _ m b, therefore NL is independent of 
0. too. 


The na —mode and the to 5—mode that are involved in the 
MI na _ m b process are characterized by a large group ve¬ 
locity mismatch \D V |, therefore the phase-matching con¬ 
dition is typically dominated by the linear term L. We 
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may therefore rewrite ^Ipk , which solves L(Qpk)+NL = 
0, as the sum of a linear contribution Qpa'.l, which solves 
L{9. pk , l ) = 0, and a nonlinear contribution f ^pk,nl << 
Qpk,l, namely LL pk =LIpk,l + Qpk,nl- 
The solution for £Ipk,l of L = 0 depends on the sign 
of D and D v . If sign(Dp3)=sign(D v )=—sign(/32) 
then we find the 2 distinct positive values flp’ K ' L . Other¬ 
wise, if sign(Dpz)=-sign(D v ) , then we find 1 positive so¬ 
lution Opjf L and if sign(D 03 )=sign(D v )=sign( 02 ) then 
no positive solutions are found. 

We can easily show that the dynamics for the MIi na - smb 
process is reversed: if sign(D p 3 )=sign{D v )=—signify) 
then no positive solutions are found, whereas if 
sign{Dp 3 )=-sign(D v ) we find 1 positive solution Qp K L 

and if sign(D^ 3 )=sign{D v )=sign{p 2 ) then we find the 2 

(1 2 ) 

distinct positive solutions U p ’ KL . 

(1 2 ) 

Therefore, overall 2 distinct solutions Ll y p ' K L are found 
for the Mlna-mb process, which can be associated to the 
MIsna-imh or the MI ina _ srnb processes. They are re¬ 


ported in Eq. (Ell¬ 
in order to calculate Qpk,nl we write f l = LIpk,l + 
Qpk,nl and we expand the powers Q 2 and f l 3 of L 
neglecting terms in Ll% K NL and &? PK NL . I n this way 
L+NL = 0 becomes a linear equation for LIpk,nl, whose 
solution is reported in Eq. By means of Qpk,nl and 
LIpk,l w e provide a fully analytical estimation of Qpk- 
Typically the bandwidth of B na - mb is small if com¬ 
pared to flpic, that is |fi — Qpk | << ^pk, and we 
can thus expand M'[l, 1] — M'[2, 2]=L(fl) + NL about 
LL pk . The derivative dL/dLl evaluated at Ll PK is 
3Dp 3 Ll 2 PK + 2/3] 2 ^pk + D v « — 2k 2 ^l~ P 1 KNL (see fourth 
line of Eq.([T7| setting LIpk,l ~ pk ), with = 
(b s /2)(C nn \p na \ 2 + C mm \p m b\ 2 ). The first order ex¬ 
pansion of L(f2) + NL is therefore L(flpp-) + NL — 
2 h 2 ^pp n UpA") — 2 /j 2 Hpp.jYp(fI and in 


Eq.(A2-l) we thus approximate (M'[l,l] — M'[2,2]) 2 = 
(L + NL) 2 « 4k%Llp 2 K nl (LI — Qpk) 2 , which results in 
the estimate of B na - mb stated in Eq. ©• 


[1] R.J. Essiambre, and R.W. Tkach, Proc. of IEEE 100, 
1035-1055 (2012). 

[2] P.P. Mitra, and J.B. Stark, Nature 411, 1027 (2001). 

[3] A.M. Rubenchik, I.S. Chekhovskoy, M.P. Fedoruk, O.V. 
Shtyrina, and S. K. Turitsyn, Opt.Lett. 40, 721 (2015). 

[4] D.J. Richardson, J.M. Fini, and L.E. Nelson, Nature 
Phot. 7, 354 (2013). 

[5] P. M. Krummrich, Opt. Express 19, 16636 (2011). 

[6] R.Ryf, R.Essiambre, J.Hoyningen-Huene, and P.Winzer, 
in Optical Fiber Communication Conference, paper 
OW1D.2 (2012). 

[7] C. Antonelli, A. Mecozzi, and M. Shtaif, in Optical Fiber 
Communication Conference, paper W3E.1 (2014). 

[8] N. Bai, E. Ip, T. Wang, and G. Li, Opt. Express 19, 
16601 (2011). 

[9] Y. Jung, S. Alam, Z. Li, A. Dhar, D. Giles, I. Giles, J. 
Sahu, F. Poletti, L. Grner-Nielsen, and D. Richardson, 
Opt. Express 19, B952 (2011). 

[10] J. Hansryd, P.A. Andrekson, M.Westlund, J.Li and 
P.O.Hedekvist,J.Sel.Top.Quant.Elec.8, 506 (2002). 

[11] R.H. Stolen, Quantum Elec. 11, 100 (1975). 

[12] A. Mussot, T. Sylvestre, L. Provino, and H. Maillotte, 
Opt.Lett. 28, 1820 (2003). 

[13] A. Tonello, S. Pitois, Wabnitz, Millot, T. Martynkien, 
W. Urbanczyk, J. Wojcik, A. Locatelli, M. Conforti, and 
C. De Angelis, Opt.Expr. 14, 397 (2006). 

[14] L. Rishj, P. Steinvurzel, Y. Chen, L. Yan, J. Demas, M. 
Grogan, T. Ellenbogen, K. Crozier, K. Rottwitt, and S. 
Ramachandran, in European Conference and Exhibition 
on Optical Communication paper Tu.3.F.2 (2012). 

[15] L.G. Wright, D.N. Christodoulides, and F.W. Wise, Na¬ 
ture Phot.(2015). 

[16] , E. Seve, P. TchofoDinda, G. Millot, M. Remoissenet, J. 
M. Bilbault, and M. Haelterman, Phys. Rev. A 54, 3519 
(1996). 

[17] S.Mumtaz, R.Essiambre, and G.P. Agrawal, J. Lightw. 
Technol. 31,398 (2013). 


[18] L. Palmieri, and A.Galtarossa, IEEE Phot. Journal 6 
(2014). 

[19] G.Millot,S.Piotis,P.T.Dinda, and M.Haelterman, IEEE 
Phot. Journal 6 (2014). 

[20] D.Modotto, G.Manili, U.Minoni, S.Wabnitz, C.De An¬ 
gelis, G.Town, A.Tonello, and V.Couderc, IEEE Phot. 
Journal 3 (2011). 

[21] G.P. Agrawal, “Nonlinear Fiber Optics (third edi¬ 
tion)”, (Academic Press,2001). 

[22] K. Nithyanandan, R.Vasantha Jayakantha Raja, K. 
Porsezian, and B. Kalithasan, Phys. Rev. A 86 (2012). 

[23] D.Anderson and M.Lisak, Opt.Lett.9,468-470 (1984). 

[24] M.Karlsson, J.Opt.Soc.Am.B12,2071-2077 (1995). 

[25] A. Labruyere, S. Ambomo, C. M. Ngabireng, P. 
TchofoDinda, K. Nakkeeran, and K. Porsezian, Opt.Lett. 
32, 1287-1289 (2007). 

[26] W.Magnus, Commun. Pure Appl. Math. 7, 649-673 
(1954). 

[27] M.Guasoni, V.V.Kozlov, and S.Wabnitz, 

J.Opt.Soc.Am.B 29, 2710-2720 (2012). 

[28] S.Buch and G.P.Agrawal, Opt. Lett. 40, 225-228 (2015). 

[29] A.Picozzi, G.Millot, and S.Wabnitz, Nature Phot. 9, 289- 
291 (2015). 




17 


[1] R.J. Essiambre, and R.W. Tkach, “Capacity Trends and 
Limits of Optical Communication Networks,” Proc. of 
IEEE 100, 1035-1055 (2012). 

[2] P.P. Mitra, and J.B. Stark, “ Nonlinear limits to the in¬ 
formation capacity of optical fibre communications,” Na¬ 
ture 411, 1027-1030 (2001). 

[3] A.M. Rubenchik, I.S. Chekhovskoy, M.P. Fedoruk, O.V. 
Shtyrina, and S. K. Turitsyn, “ Nonlinear pulse combin¬ 
ing and pulse compression in multi-core fibers,” Opt.Lett. 
40, 721 (2015). 

[4] D.J. Richardson, J.M. Fini, and L.E. Nelson, “Space- 
division multiplexing in optical fibres”, Nature Phot. 7, 
354-362 (2013). 

[5] P. M. Krummrich, “Optical amplification and optical 
filter based signal processing for cost and energy effi¬ 
cient spatial multiplexing,” Opt. Express 19, 1663616652 
( 2011 ). 

[6] R.R.yf, R.Essiambre, J.Hoyningen-Huene, and P.Winzer, 
“Analysis of Mode-Dependent Gain in Raman Amplified 
Few-Mode Fiber, in Optical Fiber Communication Con¬ 
ference, paper OW1D.2 (2012). 

[7] C. Antonelli, A. Mecozzi, and M. Shtaif, “Modeling Ra¬ 
man amplification in multimode and multicore libers,” in 
Optical Fiber Communication Conference, paper W3E.1 
(2014). 

[8] N. Bai, E. Ip, T. Wang, and G. Li, “Multimode fiber 
amplifier with tunable modal gain using a reconfigurable 
multimode pump,” Opt. Express 19, 16601-16611 (2011). 

[9] Y. Jung, S. Alam, Z. Li, A. Dhar, D. Giles, I. Giles, .1. 
Sahu, F. Poletti, L. Grner-Nielsen, and D. Richardson, 
“First demonstration and detailed characterization of a 
multimode amplifier for space division multiplexed trans¬ 
mission systems,” Opt. Express 19, B952-B957 (2011). 

[10] J. Hansryd, P.A. Andrekson, M.Westlund, J.Li and 
P.O.Hedekvist, “Fiber-based optical parametric ampli¬ 
fiers and their applications,” J.Sel.Top.Quant.Elec.8, 506 
( 2002 ). 

[11] R.H. Stolen, “Phase-matched stimulated four-photon 
mixing in silica-fiber waveguides,” J.of Quantum Elec. 
11, 100-103 (1975). 

[12] A. Mussot, T. Sylvestre, L. Provino, and H. Maillotte, 
“Generation of a broadband single-mode supercontinuum 
in a conventional dispersion-shifted fiber by use of a sub¬ 
nanosecond microchip laser,” Opt.Lett. 28, 1820-1822 
(2003). 

[13] A. Tonello, S. Pitois, Wabnitz, Millot, T. Martynkien, W. 
Urbanczyk, J. Wojcik, A. Locatelli, M. Conforti, and C. 
De Angelis, “Frequency tunable polarization and inter- 
modal modulation instability in high birefringence holey 
fiber,” Opt.Expr. 14, 397-404 (2006). 

[14] L. Rishj, P. Steinvurzel, Y. Chen, L. Yan, J. Demas, 
M. Grogan, T. Ellenbogen, K. Crozier, K. Rottwitt, 
and S. Ramachandran, “High-Energy Four-Wave Mix¬ 
ing, with Large-Mode-Area Higher-Order Modes in Op¬ 


tical Fibres,” in European Conference and Exhibition on 
Optical Communication paper Tu.3.F.2 (2012) 

[15] L.G. Wright, D.N. Christodoulides, and F.W. Wise, 
“Controllable spatiotemporal nonlinear effects in multi- 
mode fibres” Nature Phot. (2015). 

[16] , E. Seve, P. TchofoDinda, G. Millot, M. Remoissenet, J. 
M. Bilbault, and M. Haelterman, “Modulational insta¬ 
bility and critical regime in a highly birefringent fiber,” 
Phys. Rev. A 54, 3519-3534 (1996). 

[17] S.Mumtaz, R.Essiambre, and G.P. Agrawal,“Nonlinear 
Propagation in Multimode and Multicore Fibers: Gener¬ 
alization of the Manakov Equations,” .1. Lightw. Technol. 
31,398-406 (2013). 

[18] L. Palmieri, and A.Galtarossa, “Coupling Effects Among 
Degenerate Modes in Multimode Optical Fibers” IEEE 
Phot. Journal 6 (2014). 

[19] G.Millot,S.Piotis,P.T.Dinda, and M.Haelterman “Obser¬ 
vation of modulational instability induced by velocity- 
matched cross-phase modulation in normally dispersive 
bimodal fiber”, IEEE Phot. Journal 6 (2014). 

[20] D.Modotto, G.Manili, U.Minoni, S.Wabnitz, C.De Ange¬ 
lis, G.Town, A.Tonello, and V.Couderc, “Ge-doped mi- 
crostructured multicore fiber for customizable supercon¬ 
tinuum generation ” IEEE Phot. Journal 3 (2011). 

[21] G.P. Agrawal, “Nonlinear Fiber Optics (third edi¬ 
tion)”, (Academic Press,2001). 

[22] K. Nithyanandan, R.Vasantha Jayakantha Raja, K. 
Porsezian, and B. Kalithasan,“Modulational instability 
with higher-order dispersion and walk-off in Kerr media 
with cross-phase modulation,” Phys. Rev. A 86 (2012). 

[23] D.Anderson and M.Lisak,“Modulational instabil¬ 

ity of coherent optical-fiber transmission signals”, 
Opt.Lett.9,468-470 (1984). 

[24] M.Karlsson, “Modulational instability in lossy optical 
fibers”,J.Opt.Soc.Am.B12,2071-2077 (1995). 

[25] A. Labruyere, S. Ambomo, C. M. Ngabireng, P. 
TchofoDinda, K. Nakkeeran, and K. Porsezian, “Sup¬ 
pression of sideband frequency shifts in the modulational 
instability spectra of wave propagation in optical fiber 
systems”, Opt.Lett. 32, 1287-1289 (2007). 

[26] W.Magnus, “On the exponential solution of differential 
equations for a linear operator”, Commun. Pure Appl. 
Math. 7, 649-673 (1954). 

[27] M.Guasoni, V.V.Kozlov, and S.Wabnitz, “Theory of po¬ 
larization attraction in parametric amplifiers based on 
telecommunication fibers”, J.Opt.Soc.Am.B 29, 2710- 
2720 (2012). 

[28] S.Buch and G.P.Agrawal, “Soliton stability and trapping 
in multimode fibers,” Opt. Lett. 40, 225-228 (2015). 

[29] A.Picozzi, G.Millot, and S.Wabnitz, “Nonlinear optics: 
Nonlinear virtues of multimode fibre,” Nature Phot. 9, 
289-291 (2015). 



